# MEDIPS batch processing:
# Performing MEDIPS.methylationProfiling for pre-processed data

library(MEDIPS)
library(BSgenome.Hsapiens.UCSC.hg19)

getwd()
# [1] "/usr/local/data/Data/NCH/MBD-Seq/MEDIP-Analysis"

# get list of files in directory
# here we're loading the stored results of the previous batch run
files <- list.files(pattern="glioma.*RData$")

# performing the analysis for FANTOM4 promoters of chr7 
# ROI.file <- ("/home/skurscheid/Data/Riken/FANTOM4/Fantom4Promoters.chr7.txt")
# ROI.file <- ("/home/skurscheid/Data/NCH/HOX/ACGHmurat/CGH_chr7BACs.bed")

# profiling of all BACs 
# ROI.file <- ("/usr/local/data/Data/NCH/HOX/ACGHmurat/Agilent_aCGH_probes.bed.position_sorted.txt")
# output <- c("BACs.profile")

# profiling of HOXA Locus
#ROI.file <- ("/home/skurscheid/Data/UCSC/HOXA_500bp_windows.bed.txt")
#output <- c("HOXA.profile")

# profiling of CpGi in the HOXA locus
# ROI.file <- ("/home/skurscheid/Data/UCSC/HOXA_CpGi.txt")
# output <- c("HOXA.CpGi.profile")

# profiling of the whole genome
chrom <- "chr"

#####################################
# beginning of loop
#####################################

#loop through files, process and write results to individual RData files 
for (i in files) {
	sample <- gsub(".RData","",i)
	print(sample)
	
	# loading the RData file
	load(i)
	
	# using median of read counts for the methyProfiling function to better compensate for outliers.
	# here for a ROI file	
	# profile <- MEDIPS.methylProfiling(data1 = MEDIPS.SET, ROI_file = ROI.file, math = median, select = 2)
	
	# loop for processing all autosomes
	for (j in (1:22)) {
	  chrom <- paste("chr",j,sep="")
	  print(chrom)
	  output <- i
	  profile <- MEDIPS.methylProfiling(data1 = MEDIPS.SET, chr=chrom, frame_size=500, step=250, math = median, select = 2)
	 	# extraction of AMS (Absolute Methylation Score) values, as these are corrected for CpG densities (Coupling Factors)
	  # extraction of RPM (Reads Per Million)
	  # extraction of RMS (scaled Reads Per Million - log2 transformed [0:1000])
	  if (which(chrom == j) == 1){

	  } else {
		  profile.ams <- cbind(profile.ams, as.matrix(profile$ams_A))
		  colnames(profile.ams)[which(chrom == j)] <- chrom
		  profile.rms <- cbind(profile.rms, as.matrix(profile$rms_A))
		  colnames(profile.rms)[which(chrom == j)] <- chrom
		  profile.rpm <- cbind(profile.rpm, as.matrix(profile$rpm_A))
		  colnames(profile.rpm)[which(chrom == j)] <- chrom
	  }
	  # write collated profiling results to CSV
    write.csv(profile.ams, file=paste(output,"ams.csv", sep="."))
    write.csv(profile.rms, file=paste(output,"rms.csv", sep="."))
    write.csv(profile.rpm, file=paste(output,"rpm.csv", sep="."))
	}
	gc()
}
	
## write collated profiling results to CSV
#write.csv(profile.ams, file=paste(output,"ams.csv", sep="."))
#write.csv(profile.rms, file=paste(output,"rms.csv", sep="."))
#write.csv(profile.rpm, file=paste(output,"rpm.csv", sep="."))

###########################################
# end of loop
###########################################

